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Abstract: We propose a new technique to measure the time delay of radio- loud gravitational 
lens systems, which does not rely on the excessive use of interferometric observations. Instead, 
the method is based on single-dish flux density monitoring of the (unresolved) lens system's total 
lightcurve, combined with additional interferometric measurements of the flux density ratio at 
a few epochs during that monitoring period. 

The basic idea of the method is to reconstruct the individual image lightcurves from the 
observed total lightcurve by assuming a range of potential values for the time delay and the 
magnification ratio of the images. It is then possible to single out the correct reconstruction, and 
therefore determine the time delay, by checking the consistency of the reconstructed individual 
lightcurves with the additional interferometric observations. We performed extensive numerical 
simulations of synthetic lightcurves to investigate the dependence of the performance of this 
method on various parameters which are involved in the problem. Probably the most promising 
candidates for applying the method (and also for determining the Hubble constant) are lens 
systems consisting of multiply imaged compact sources and an Einstein ring, such as B0218-I-357 
from which some of the parameters used for our simulations were adopted. 

1 Introduction 

Already many years ago it was realised (Refsdal 1964) that gravitational lensing provides a way 
to determine the Hubble constant. This method is based on measuring the time delay At between 
the arrival times of light rays corresponding to different images in gravitational lens systems. 
Of course, an intrinsic variation of the source on appropriate timescales is essential for the time 
delay to be detectable. In addition, the redshifts of source and deflector and a reliable model 
for the mass distribution of the deflector are required to calculate the Hubble constant. The 
advantage of the time delay method (as also of the Sunyaev-Zel'dovich effect and potentially also 
of SN la) compared to conventional techniques for determining Hq is the possibility to reach 
cosmological distances directly without using a "distance ladder" of intermediate calibration 
steps. 

The first lens system discovered, the double quasar 0957-1-561 (Walsh et al. 1979), has 
been monitored since then at optical and radio frequencies, and quite a number of papers 
about its time delay have been published. However, the value of At still remains controversial 
today, because these studies were plagued by problems like unevenly sampled data points and 
microlensing effects. For an overview of the status quo see Chap. 2 of Kochanek & Hewitt (1996). 
But even if a secure value for At were available, this would not be of much use for determining 
the Hubble constant, because it is extremely difficult to specify a unique model for the lens mass 
distribution in this system. This problem arises because not only a single galaxy, but also the 
cluster of galaxies in which it is located, as well as another cluster at a different redshift, are 
contributing to the image splitting, and so even in the simplest realistic models there are more 



parameters than can be constrained observationally (e.g. Kochanek 1991, Bernstein et al. 1993; 
but also see Grogin & Narayan 1996). 

A very much better candidate for determining Hq is the lens system B0218+357 (Patnaik et 
al. 1993), consisting of two images of a compact flat spectrum radio source and an Einstein ring 
of extended emission. This system combines several important advantages. The small image 
separation of 0".335 implies that probably a single galaxy is acting as a lens. In addition, the 
morphology of the compact images (Patnaik et al. 1995) and especially the radio ring provide 
valuable information about the mass distribution of the lensing galaxy. Thus it should be feasible 
to construct a sufficiently accurate and well-constrained lens model for this system. Furthermore, 
as the time delay for this system is expected to be of the order of two weeks, and the source 
is known to be variable, a determination of At should be possible on reasonable timescales of 
a few months. In fact, Corbett et al. (1996) already published a value of At = 12 ± 3 days, 
derived from the time variation of the image polarization in VLA monitoring observations. 

In this paper we are investigating the possibilities for determining the time delay At of 
radio-loud gravitational lens systems without the excessive use of expensive interferometric ob- 
servations. Instead, we consider monitoring the total flux density of the unresolved lens system 
with a single-dish radio telescope. As we briefly motivate in Sect. ^ an analysis of the autocor- 
relation function of the combined lightcurve is a straightforward approach to extract the value 
of the time delay. However, numerical simulations show that it is difficult to apply this method 
in practice, because it would require unrealistically long monitoring periods in order to achieve 
significant results. Therefore the main part of this paper (Sect. ^) is devoted to a new method 
which is shown to yield a much more reliable result for the time delay by making use of a few ad- 
ditional interferometric measurements. The idea of this method is to reconstruct the lightcurves 
of the individual images from the observed combined lightcurve by assuming values for the time 
delay and the magnification ratio of the images. To single out the true value of the time delay we 
use the consistency of the reconstructed lightcurves with additional information about the flux 
density ratio of the images at different epochs obtained from interferometric observations. In 
extensive simulations we investigate the performance of this "lightcurve reconstruction method" , 
depending on the various parameters involved, and show that typically a handful of interfero- 
metric observations suffice to determine At reliably. Finally, in Sect. Qwe summarize the results 
and discuss the prospects for the application of this new method, especially in context with the 
most promising lens system B0218-I-357. 

The discussion in this paper is restricted to gravitational lens systems with two images of a 
compact source, although the formalism described here might as well be extended to a larger 
number of images. As we are mainly interested in radio lightcurves it is justified to neglect 
microlensing effects for this study. Furthermore we assume that there is no variation on relevant 
timescales in the flux density of extended emission which may be associated with the compact 
source (e.g. the ring in B0218-I-357). 
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2 Autocorrelation approach 



2.1 Time delay effect on the autocorrelation function 

The total lightcurve S{t) of a two-image gravitational lens system results from a time shifted 
superposition of the intrinsic source lightcurve Sin{t) and a constant contribution ^const (e.g., 
from an extended source component) in the following way: 

S{t) = Si{t) + S2{t) + 5const = AilSin(t) + l^2Sin{t - At) + 5const . (1) 

5i_2(t) and fii^2 denote the flux density and the absolute magnification of the individual images.[] 
Introducing the observable magnification ratio fi := /^i//^2 and subtracting the time averages of 
the flux densities, S{t) := S{t) - {S{t))t, Si{t) := Si{t) - {Si{t))t and S2{t) := S2{t) - {S2(,t))t, 
Eq. (H) becomes 

S{t) = Slit) + S2{t) = Slit) + -Slit - At) . (2) 

A* 

Here we assume that 5(t) is a stationary random process and deflne its (normalized) autocorre- 
lation function C(r) as 

Cir) := III with e(r) := (5(t) 5(t + r))^ . 

It is obvious that the time delay should show up in C(r) as a positive contribution to the 
autocorrelation at r = At. Starting from Eq. (^) it is an easy exercise to show that C(r) can 
be calculated from the intrinsic autocorrelation function Ci^ir) of the source according to 



Cir) = - 
a 



Cm(^) + ^^C'in(r - At) + -Y^Q^ir + At) 
+ 1 /i"^ + 1 



(3) 



with the factor a = 1 + ^^^-^ Cin(At) ensuring the proper normalization. Corresponding to the 
change in the autocorrelation function the power spectrum is modified according to 



P(u;) =^?Pin(w) 



1 + ^ + - cos(wAt) 



(4) 



where Pini^^) denotes the intrinsic power spectrum of the unlensed source. 

Figure |l|a graphically depicts the modification of the autocorrelation function described by 
Eq. dll). Here we took the values At = 20 (in arbitrary time unitsQ) for the time delay, ^ = 3 
for the magnification ratio and 

2 2 

CUr) = -var^;^^^ (5) 

as an example for the intrinsic autocorrelation function, which corresponds to the reasonable 
intrinsic power spectrum 

Pin(w) oc we""^™^ . (6) 

The parameter rvar represents a characteristic timescale of the intrinsic source variability and has 
been chosen to be Tyar = 3.18 for this plot. The resulting autocorrelation function C(t) shows 

^In the following we will restrict the time delay to be positive. Hence the index 1 denotes the image in which 
intrinsic source variations will appear first. 

■^By scaling all time variables appropriately the arbitrary time units used here can be adapted to any lens 
system of interest. 
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Figure 1: (a) The autocorrelation function C(r) (solid 
line) resulting from a time shifted superposition of an 
intrinsic lightcurve with autocorrelation function dnir) 
(dashed line). The value of At is indicated by the dot- 
ted line, (b) This diagram shows the autocorrelation 
function calculated for a simulated lightcurve with pa- 
rameters as shown in Table |^. Again the dotted line 
indicates the value for the time delay. 

a distinct maximum at the value of the time delay. Thus it is at least in principle possible to 
obtain information about At by analysing the autocorrelation function of combined lightcurves. 
However, if the timescale of variability Tvar is comparable to or longer than At, the "time delay 
peak" will merge with the "intrinsic maximum" of C(t) at t = and any information on At 
will be lost. 

2.2 Analysis of C(r) for simulated lightcurves 

Of course, in applications to observations the limited observing time T, the finite sampling 
interval AT and observational errors will impose serious constraints on the usefulness of the 
autocorrelation function for determining the time delay. To study this quantitatively we used 
synthetic data sets of combined lightcurves. These lightcurves were generated as realizations 
of a gaussian random process with the intrinsic power spectrum of Eq. (^). For simplicity we 
restricted the study to constant observing intervals AT. 

In Fig. I^b the autocorrelation function for one typical lightcurve realization with parameters 
as summarized in Table |^ is plotted. Due to the limited observing period T, which does not 
provide a "fair sample" of the lightcurve's statistical properties, a number of additional maxima 
and minima are now showing up, making it very difficult to identify the time delay peak. A 
simple quantitative measure for the performance of the autocorrelation method is the fraction 
V of lightcurve realizations for which the time delay maximum, i.e. the maximum closest to the 

Table 1: Standard parameter values for simulated lightcurves. The parameter rj is defined as the ratio of the 
standard deviation ass of observational errors added to every data point and the dispersion as of the lightcurve 
itself. 



parameter 
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standard value 


observing period 


T 


100 


observing interval 


AT 
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V 


0.1 


power spectrum 
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timescale of variability 
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Figure 2: The dependence of V (see text for definition) 
on (a) the observing period T and (b) the magnification 
ratio ^ of the images; all other parameters are as in 
Table 0. 



true time delay, is in fact the highest maximum within the range of potential At values. In the 
following all maxima in the range [0, 50] except the r=0 maximum are considered. We studied 
the dependence of V (and other measures) on all of the parameters shown in Table |l|, but here 
we only present the results for varying T and ^, keeping the other parameters fixed. 

From Fig. ^a it can be seen that prolonging the observation period T leads to an increase 
in because with increasing T the sampling of the lightcurve improves and the observed 
autocorrelation function converges to the theoretical C(r) shown in Fig. |^a. However, since 
for some realizations the unwanted additional maxima in C(r) can be quite stable, it takes a 
rather long observing time to ensure getting significant results for the time delay. In Fig. ||b 
the dependence of V on the magnification ratio is depicted. Of course the method works best 
for = 1, as for increasing magnification ratio the combined lightcurve is dominated more 
and more by the brighter image. The autocorrelation function C(r) will then be dominated by 
the first term in Eq. (|^), which does not contain the time delay, and converges to the intrinsic 
autocorrelation function Cin(T) for /i oo. 

To summarize, one can say that an analysis of the autocorrelation function of the total 
lightcurve of gravitational lens systems only constitutes a viable method for determining the 
time delay if the magnification ratio is close to unity and the observing period T is very much 
longer than At and Tvar- 



3 Lightcurve reconstruction method 



In this section we introduce the "lightcurve reconstruction method" for measuring the time 



delay of gravitational lens systems. As we will explain in Sect. |3.1| , it is possible to reconstruct 
the (unobserved) lightcurves Si{t) and 52 (t) of the individual images, if we assume values for 



the time delay At and the magnification ratio /U. Sect. |3.2| shows qualitatively how the true At 
value can then be singled out by checking the consistency of the reconstructed lightcurves with 
additional information from interferometric measurements about the flux density ratio of the 
images. In Sect. 3^ this is discussed quantitatively in terms of the function. Using lightcurve 
simulations we investigate in Sect. ^.4| the dependence of the results on the various parameters 



that are involved in this problem. Finally, Sect. 3.5 shows how confidence intervals for the time 
delay can be obtained in individual realizations. 
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Figure 3: The principle of lightcurve reconstruction: (a) The observed combined hghtcurve. (b) Reconstruction 
of the individual lightcurves with Atji = 20 and /xr = 2 by iteratively applying Eq. ^ (solid arrows) and 
Eq. ^ (dashed arrows), (c) The errors cr[S']^j^] and o"[S2r] of the reconstructed individual lightcurves under the 
hypothesis that Atn and /xr are in fact the true values. These were calculated from Eqs. ( |l^ and ( p^ ) with 
constant observing error (t[5^"'] = ass = ??o"s and using an estimate for the initial errors as described in the text. 



3.1 Reconstruction of lightcurves 

It is convenient to discuss the method by looking at Fig. ^. As an example for an "observed" 
combined lightcurve Fig. 0a displays one realization of a gaussian random process with param- 
eters as shown in Table |l] except that here we took fi = 2. Now, to reconstruct the individual 
image lightcurves S'iR(t) and S'2R(t) we have to postulate values At^ and fj,^ for the time delay 
and the magnification ratio, respectively.^ It will become clear in the following subsections how 
to actually determine these quantities with this method. In the example of Fig. |3| we used just 
for illustration the values Atpj, = 20 and //r = 2 which in fact agree with the "true" values used 
for generating the combined lightcurve. In addition, we have to make a guess about the individ- 
ual lightcurve S2R{t) in the time interval t € [0, AtR[, but luckily the choice for this initial guess 
will turn out to be rather irrelevant for the method to work. For simplicity we use S2R{t) = 
for t G [0, AtR[ within this paper. 

Having specified Atp,, /iR and S2R{t) for t € [0, AtR[ we can immediately calculate the 
individual lightcurve Sinit) in the same interval using 

Smit) = S{t) - S2Rit) (7) 

which is trivially derived from Eq. (^). The application of Eq. (|^) is visualized by a solid arrow 
■^Here and in the following the subscript r indicates quantities used for or obtained from the reconstruction. 
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in Fig. ^b. The next step is to compute 52R(t) on the subsequent interval [A%, 2AfR[ using 

S2R{t) = —Sinit - Atn) • (8) 

This is indicated by the dashed arrows of Fig. ^3. By iteratively applying Eqs. ^ and the 
individual image lightcurves can be reconstructed for the total observing period. Introducing 
the notation S^"^ for the restriction of the function S{t) to the interval [{n — 1) Atpj,, n AtR[, the 
reconstruction process can be written asQ 

c,{n) _ I Q(n-m) , I '-> qW (q\ 
^2R - 2^ , n-1 •^2R 



m=l 

and 



/^R /«R 



„{n) _ sr^ Q{n-m) , [-^) n{l) (-> p.\ 

- 2^ ,,m + ,,n-l ' i^^Uj 

m=0 /^R 

Now we would like to quantify the rms errors and (^[S'jp^] of the reconstructed indi- 

vidual lightcurves under the hypothesis that the values At^ and used for the reconstruction 
are in fact the true ones. Therefore these errors include the observational errors of the total 
lightcurve measurements and the errors caused by the arbitrary initial guess. In the following 
we assume the observational errors to be normally distributed with standard deviation cj[5''-"''] 
and statistically independent for each data point. Hence the increase of with respect to 

'^['5'2r] when applying Eq. (0) is described by 

a2[5S]=a2[5(^)] + a2[5(")], (11) 
whereas using Eq. (^) leads to an evolution of the error according to 

a[s'^^'\ = -a[S^;S]. (12) 
/^R 

In analogy to Eqs. (^) and (0) the error propagation during the reconstruction can be compactly 
expressed as 

n-l -. 1 

-^[^R^] = E ^LAs^--^ + ^As^i] (13) 

=1 ^R AiR 



m= 



and 



AS'^\ = E + -^,AS^i] ■ (14) 

m=0 /"r 

However, looking at Fig. |3|c is much more illustrative to understand the evolution of the errors 
of the reconstructed lightcurves. For the first interval [0, AtR[ we have to estimate the error 
C7[S2^] caused by the arbitrary initial guess. Here we assume that the variability characteristics 
of the second lightcurve in the starting interval does not substantially differ from the observed 
variability of the combined lightcurve during the total observing period. This leads to the 
estimate o"^[S'2r] ~ "rr~2"'^i '^ith ^5 denoting the dispersion of the combined lightcurve during 



■'Note that the quantities in this notation are still meant to be functions of t although the argument has been 
omitted. 
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the observation period.^ Figure ^ clearly shows that the initial uncertainty in the reconstructed 
lightcurves decreases at every interval, because the errors are multiplied with the factor 1//Upj 
in Eq. (|l2[). Thus the reconstruction is most accurate at the end of the observing period. Here 
the errors are dominated by the observational errors of the combined lightcurve, which are 
added when applying Eq. (11), and the initial uncertainty has almost dropped out. Of course 
the reconstruction works best for large values of the magnification ratio /iR because then the 
initial errors are decreasing faster. In the following we will restrict //r to be larger than one, 
because otherwise the errors will actually grow and the reconstruction fails. However, this is no 
fundamental problem for the method, since for fi^i < 1, i-e. in cases in which the light arrives 
first in the weaker image, the reconstruction can be done starting at the end of the observing 
period. 

For simplicity we will assume for this study that the standard deviation of the observational 
errors of the combined lightcurve is constant for all data points, (t[S'*^"'^] = a^s- In addition, we 
restrict the simulations performed in this paper to constant observing intervals AT. It should 
be emphasized here that the limitation to constant ass is merely for convenience, and 

removing these restrictions does not impose any conceptual difficulties for the reconstruction 
method. An important point to be mentioned is the following. The errors c'"[S'ipj,] and o"[5'2r] 
shown in Fig. |^ are strictly correct only for integer multiples of the observing interval AT. 
In order to calculate the reconstructed lightcurves for t values lying between the data points, 
the observed combined lightcurve has to be interpolated which in general introduces additional 
errors. Interpolation also is necessary to do the reconstruction with Atji not being an integer 
multiple of AT. We use a linear interpolation and determine the error introduced by the 
interpolation directly from an ensemble of simulated lightcurves. For real data, however, an 
extrapolation of the variability characteristics to timescales shorter then the sampling interval 
AT has to be made in order to get an estimate of the errors introduced by interpolation. In 
practice this could be achieved by following the approach of Press et al. (1992a). They describe 
a method that provides a "reconstruction of a set of irregularly sampled measurements into 
a continuous function and an associated standard error function" by using an estimate of the 
underlying autocorrelation function obtained from the data. 

Figure Q shows another example for a lightcurve generated according to a gaussian random 
process^ with the parameters of Table ||. In Fig. ^ the unobserved individual image lightcurves 
Si{t) and 5*2 (i) are depicted. With At = 20 and fi = 3 these add up to the observed combined 
lightcurve S{t) in Fig. In Fig. ^ the individual lightcurves are reconstructed using the true 
values for AtR and ^r, whereas Fig. ^ shows the reconstruction using the wrong value AiR = 10 
for the time delay. In agreement with the error discussion above the correctly reconstructed 
lightcurves resemble the true ones quite closely towards the end of the observing period. Of 
course now the problem arises how to distinguish between the correct reconstruction (Fig. ^c) 
and the wrong one in Fig. ^ as well as all the other wrong reconstructions conceivable with 
different AtR's and ^r's. Although it might be possible in some special cases to infer At by 
making assumptions about the shape of correctly reconstructed lightcurves, this is certainly not 
viable in general. Therefore we clearly need some additional information to solve the problem. 



^For lightcurves with timescale of variability comparable to (or larger than) Ata the dispersion is increasing 
when extending the interval on which it is calculated from AtR to T. Therefore, the dispersion used for this 
estimate should be calculated on an interval equal to AtR in such cases, because otherwise the error of the initial 
guess would be overestimated. 

®Note that the reconstruction method does not rely on properties like gaussianity or stationarity. The reason 
to use gaussian random processes for the simulations is just because they are easy to generate. 
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Figure 4: One realization of a gaussian 
random process generated with param- 
eters according to Table 0. (a) The 
unobserved individual lightcurves Si{t) 
and S2{t). (b) The observed combined 
lightcurve that results by adding up Si (t) 
and 52 (t) with At = 20 and fj. = 3. (c) 
The reconstructed individual lightcurves 
Smit) and S2ii{t) using the true values 
AtR = 20 and fiR = 3. (d) The re- 
construction when using fiR = 3 and the 
wrong value AtR = 10 for the time delay. 



3.2 Flux density ratios as additional constraints 

For the rest of this paper we investigate the chances for determining the time delay by using 
additional observational constraints about the flux density ratio of the images obtained from 
interferometric observations. It is useful to take the flux density ratio rather than absolute flux 
densities of the images, because measuring the former does not suffer from calibration problems 
and is therefore much more accurate. 

The time dependent flux density ratio m{t) of the two images is given by 

^ {S,{t)), + S^{t) ^ (^i+2)T^ + 5i(t) 

S2{t) {S^m + S2{t) (5i+2)T^ + 52(t) ' 

where Si{t) and 5*2 (t) denote the flux densities without subtraction of the mean values (cf. the 
definition in Sect. In the final identity of this equation the individual image time averages 
have been expressed in terms of the combined lightcurve average (5*1+2) := {Si{t) + S2{t))t- Note 
that (5*1+2) only includes the flux densities of the varying components in the lens system. Thus 
in order to determine (5*1+2) at least one interferometric observation is required to subtract the 
constant flux density of non varying components from the average single-dish flux density. 

In analogy to Eq. ( [T5| ) we can calculate the reconstructed flux density ratio m^{t) from the 
reconstructed individual image lightcurves as follows: 

(^1+2) + SMt) (^1+2) + S{t) - S2K{t) 
mnit) = = = ; . (16) 

(5l+2)i4^ + 52R(t) (5i+2)^ + 52R(t) 
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In order to express TnYi{t) in terms of quantities with independent errors we replaced Siji{t) with 
S{t) — S2R{t) in this equation, because for a given epoch ti the reconstruction error cr[S2R(ti)] 
only contains observational errors of S{t) for epochs preceding tj and therefore is independent 
of a[S{ti)]. By applying the usual error propagation law we obtain the relation 



CT2[mR(t)] 



52R(t) 



1 



Sm{t) 



(5*1+2) 147^ + S2R{t) 



{S1+2) + Smit) + S2n{t) 



{Sl+2)i^ + S2R{t) 



(17) 



for the error cr[mR(f)] of the reconstructed flux density ratio under the hypothesis that AtR and 
yUR are the true values. 

Figure ^ shows the unobserved true time evolution of the flux density ratio m{t) for the 
example of Fig. ^ as well as its correct and wrong reconstructions mR(t), corresponding to the 
lightcurve reconstructions in Fig. ^ and d. In order to compute the flux density ratio curves for 
this example we have to specify the value of (;Si+2) appearing in Eqs. (|^) and (|l^). To do so 
we introduce a new parameter for our lightcurve simulations, which is not included in Table |l] 
because it was not relevant for the autocorrelation analysis. We define the parameter 



(^1+2) 



as the ratio of the total lightcurve's dispersion to the average combined flux density in order 
to describe the relative variability of the source. For the example case a variability of = 2% 
has been chosen. The reconstructed flux density ratio curves in Fig. ^ also include la error 
bars calculated from Eq. (p!7|). Again it can be seen that the initial uncertainty gradually 
decreases during the reconstruction process. At the end of the observing period (t = 100) the 
contributions of ^"[^(i)] and o"[S'2r] to the error of the flux density ratio are comparable, whereas 
that of a[{Si+2)] is negligible. 

From Fig. ^ it is clear that an additional measurement of the flux density ratio, say at 
ti = 100, can rule out the combination of the values A^r = 10 for the time delay and /xr = 3 



(a) 



(c) 





100 



Figure 5: (a) The unobserved time evo- 
lution of the flux density ratio m{t) for 
the example lightcurve plotted in Fig. ^ 
(b) The flux density ratio m^,{t) recon- 
structed with the true values At^ = 20 
and /iR — 3. (c) The wrong reconstruc- 
tion with At-R = 10 and /ir — 3. 
The filled circles ai t = 100 and the dot- 
ted vertical lines in the top panel indicate 
additional interferometric flux density ra- 
tio measurements. 
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Figure 6: (a) The hatched area includes all pairs of values for Atn and fin with the reconstructed flux density 
ratio mj^(ti = 100) being compatible with the actually observed value m{ti = 100) within a 3a error interval. In 
(b) and (c) the regions of 3a consistency with flux density ratio constraints at t2 = 95 and tg = 90 have been 
added. A^i denotes the number of interferometric measurements. Corresponding to the diagram in (c), (d) shows 
the 3a contours of the function. 



for the magnification ratio, because these would lead to a completely wrong prediction for the 
flux density ratio at that epoch. However, we have to test the compatibility of all potential 
values for AtR and /xr with the additional constraint m{ti = 100). This has been done in 
Fig. ^a. The hatched area in this plot includes all pairs of parameter values (AtR, //r) with the 
reconstructed flux density ratio 77iR(ti = 100) being compatible with the measured flux density 
ratio m{ti = 100) within a 3a error margin, i.e. 

|mR(ti) - m{ti)\ < 3 ^ a'^irriRiti)] + a^[m{ti)] . 

The values for the flux density ratio measurements m{ti) used here and below include a normally 
distributed observational error with standard deviation a[m{ti)] = 0.005 m(tj). Note that for 
any time delay AtR there is, at least for this example, a value for the magnification ratio /^r 
which yields a reconstruction in agreement with the additional flux density ratio constraint. 
Thus in order to restrict the range of possible time delay values, more than one flux density 
ratio measurement is required. Figure ^3 additionally shows the 3cj consistency region defined by 
another flux density ratio measurement performed at t2 = 95, and in Fig. |6|c a third measurement 
at ta = 90 is included. The intersection of these regions is fllled black. These plots nicely 
illustrate how the method works. However, to quantify the consistency of parameter values 
(AtR, /iR) with the interferometric constraints a discussion in terms of is adequate. 

3.3 Analysis of minima in the function 

In general functions are used to flt parametrized models to a given data set and to quantify the 
compatibility of the model with the data. In the case discussed here the "model" is represented 
by the basic assumptions made for this study which were formulated in the last paragraph of 
the introduction. The postulated values AtR for the time delay and /xr for the magnification 
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ratio can be regarded as the fit parameters. Therefore we define as a function of Atj^ and 

2/A+ ^ ^ [mRitj ; A^r, /in) - m{ti)f 



A^i denotes the number of interferometric flux density ratio measurements. Note that here not 
only the observational data m(tj) but also the values mR(tj ; AtR,//R) predicted by the model 
contain errors^ which have to be included in the x^ function. Another peculiarity of the x^ 
used here is the extreme non-linearity of the model in its fit parameter AtR. This leads to the 
ocurrence of a variety of minima in the x^ function and to non-elliptic confidence regions for 
the parameter values. However, as we will see below the model is still sufficiently linear locally 
(at individual minima), so that x^ statistics can be applied (see e.g. Sect. 15.6 of Press et al. 
1992b). 

Figure ^ shows the 3a contours of the x'^ function for the same number of interferometric 
measurements as in Fig. ^c. The definition of the 3a consistency value for x^ which has been 
used here will be given below. This plot reveals that the information from the three flux density 
ratio measurements included here does not suffice to determine the time delay unambiguously. 
There are still a number of 3a minima left. In the following we will call the minimum which 
is closest to the true parameter values the true minimum and all other minima will be termed 
wrong minima. Figure |^ depicts the evolution of the 3a contours when the number A'^, of 
interferometric measurements is further increased. It can be seen here that the wrong minima 
are gradually vanishing. For Ni = 6 the true minimum is the only minimum remaining and 
therefore a definite time delay determination can be obtained. However, already for Ni = 4 the 
global minimum in the x^ function, which gives the best fit parameter values, corresponds to 
the true minimum. 

To test the lightcurve reconstruction method we again used simulated lightcurves with pa- 
rameters as shown in Table |l|. Additional parameters which have to be specified for the lightcurve 
reconstruction are displayed in Table This table also shows the "standard values" that we 
have adopted for these parameters in the simulations. For simplicity the intervals between the 
interferometric observations are assumed to be constant. These observations are performed to- 
wards the end of the observing period T when the reconstruction is most accurate. Again we 



^These errors also depend on the fit parameters Atn and /ir. Strictly speaking, this parameter dependence 
of the errors leads to an additional term which has to be considered when deriving the function from the 
maximum likelihood principle. However, this term can be neglected here because the variation of the errors is 
small for parameter values of interest. 



3.2 r 

1 3.0 r 

2.8 r 
2.6 L_ 

3.2 r 
a" 3.0 r 
2.8 r 
2.6 L. 

3.2 r 
:l 3.0 r 
2.8 r 
2.6 E_ 



10 



20 



30 



N,=4 



N,=5 



N,=6 



40 



50 



Figure 7: These plots show the evolution of the Scr con- 
tours when increasing the number Ni of interferometric 
observations. For A^i = 6 only the true minimum re- 
mains. The global minimum of the function within 
the parameter range considered is marked with a cross. 
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Table 2: Standard parameter values used for applying the lightcurve reconstruction method to the simulations. 



parameter 


symbol 


standard value 




7/ 


2% 


Timnhpr of iiitprfprnTTiptTir' ohsprvi^tioTis 




5 


interval between interferometric observations 


ATi 


5 


error of flux density ratio measurements 


a[m]/m 


0.5% 


error of the average combined flux density 


a[(5i+2>]/(5i+2) 


2% 



considered all minima in the range [0, 50] as potential time delay values.^ Figure ^ summarizes 
the results of applying the reconstruction to 1000 simulated lightcurves with the parameter val- 
ues of Tables || and Figure ^ shows the cumulative distribution of the AtR values at the 
global minimum of the function. From this plot it can be seen that for most of the simulated 
lightcurves the location of the global minimum in fact agrees with the value for the time delay 
which has been used for generating them. To get an estimate of the accuracy achievable for the 
time delay determination, Fig. ^ depicts the distribution of the At^i values at the true mini- 
mum. The narrow distribution centered on the exact value of At = 20 shows that a fairly high 
precision can be expected with the standard parameter values used here. However, to obtain 
error ranges for individual realizations, statistics can be applied. This will be discussed in 
Sect. |3^. Plots c and d of Fig. |8| show the distribution of the /iR values at the global and at 
the true minimum, respectively. For the rather small lightcurve variability of u = 2% it is not 
surprising that the magnification ratio is well determined here. The two distributions are nearly 
identical, because also in the few cases where the global minimum is a wrong minimum, it is 
located roughly at the correct value for the magnification ratio (cf. Fig. ^d). Figure ^ shows 
the distribution of the rank of the true minimum when ordering the minima according to rising 
values. For 93% of the lightcurves the global minimum in the function corresponds to 

*In practice, of course, a plausible range of potential time delay values can be given if a lens model is available. 
For B0218+357, e.g., the time delay is expected to be roughly between 8 and 20 days. 
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Figure 8: Plot (a) shows the cumulative 
distribution of At-R, values at the global 
minimum of the function, whereas in 
(b) the distribution of At^ values at the 
true minimum is depicted. Accordingly 
plots (c) and (d) show the distribution of 
the fiR values at the global and at the true 
minimum. Diagram (e) is a histogram of 
the rank of the true minimum when or- 
dering the minima according to rising 
values. Finally, (f) shows the distribu- 
tion of the number of Scr minima appear- 
ing in the x^ function. 
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the true minimum, and only for ~ 1% the rank of the true minimum is higher than 2. Finally, 
Fig. 1^ displays the distribution of the number of 3cr minima which are appearing in the ^ 
function. In the majority of cases just one or a few minima are occurring, but there are also 
rare cases with a large number of minima showing up. 

In the following paragraph we will check if some statements given by ^ statistics are ap- 
plicable for this analysis. Let denote the number of constraints used in a fit and n be 
the number of fit parameters. According to statistics the probability density distribution for 
the value at the minimum of the x^ function is given by a -distribution with m := N — n 
degrees of freedom, i.e. 

Px^mix ) = Y{^m)2T ^ ^ e ^ . 

This is valid if the errors are normally distributed and the model is linear in its fit parameters. 
In our case, however, the model is extremely non-linear in the fit parameter At^, and as we 
have mentioned before this leads to the occurrence of additional wrong minima which cannot 
be described by x^ statistics. As we have seen in Fig. ^e the wrong minima can be lower 
than the true minima. Therefore the relation given above is certainly not valid any more for 
the global minimum of the x^ function. Nevertheless the x^ values at the true minima still 
obey a x^-distribution with m degrees of freedom. This is visualized by the thick, solid-line 
histogram in Fig. ^ which shows the distribution of the values at the true minimum for the 
simulated lightcurves and which matches quite well with the overlaid theoretical x^-distribution 
with 5 — 2 = 3 degrees of freedom as it is expected for A'l = 5 constraints and 2 fit parameters. In 
addition, this figure reveals that the distribution of the x^ values at the exact parameter values 
is consistent with a x'^-distribution with A'l = 5 degrees of freedom, again in agreement with the 
predictions of x^ statistics. The thin, solid-line histogram in Fig. ^ depicts the distribution of 
the x^ values at the lowest wrong minimum of the x^ function which of course is not described 
by x^ statistics. The dotted vertical lines in the figure indicate limiting x^ confidence values 
for rejecting the hypothesis that some given parameter values AtR and /^r are the correct ones. 
These exclusion confidence values have been derived from the x^"distribution with A'^i = 5 
degrees of freedom for x^ at the exact parameter values. Corresponding definitions already have 
been used for the So" contours in Figs. ^ and and for quantifying the frequency of wrong 
minima in Fig. ^f. Actually this is a conservative limit for rejecting wrong minima, because for 
the exclusion of the hypothesis that a given minimum is the true one we would have to use the 

\a 2a 3ct 



Figure 9: The thick, sohd-line histogram shows the dis- 
tribution of values at the true minimum overlaid 
with the theoretical probability density distribution 
with 5—2 = 3 degrees of freedom. The dotted histogram 
represents the distribution of values at the exact pa- 
rameter values At = 20 and /i = 3 overlaid with the 
theoretical distribution with 5 degrees of freedom. 
The dotted vertical lines indicate the x^ limits for a 
given exclusion confidence according to the theoretical 
distribution with 5 degrees of freedom. The thin, solid- 
line histogram displays the distribution of values at 
the lowest wrong minimum. (This distribution is ex- 
tending beyond the range of x^ values covered by this 
plot). 
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X^-distribution with Ni — 2 degrees of freedom which is shifted to smaUer values. In Sect. 3.5 
we will show with an example how to assign significance values to the statement that one of 
several minima in the function is indeed the true minimum. 



3.4 Parameter dependence of the results 

In order to investigate the parameter dependence of the performance of the lightcurve recon- 
struction method we introduce the following measures to quantify the information which for the 
standard parameter values is contained in Fig. ^: 

• The probability "Pj, i = 1, 2, 3 for the true minimum being the global, the second lowest or 
the third lowest minimum in the function (cf. Fig. ^e). 

• The average number Mj, i = 1, 2, 3 of lex, 2a and 2>a minima appearing in the x^ function 
(cf. Fig. |f). 

• Estimates a{/S£) and (t(^) for the errors of the position of the true minimum derived from 
the distributions of the AtR and /xr values at this minimum (cf. Fig. |b and d). 

In the following we study the dependence of these quantities on the parameters shown in Tables |l] 
and 0. Each time one parameter is varied and the others are kept fixed at the standard values. 



Time delay At: Figure IC presents the results for changing the value of the time delay At 
between and 50. This corresponds to the Atpj range which is considered as the range of 
potential time delay values and for which the x^ function is calculated. The graphs indicate 
that apart from some boundary effects the results are rather insensitive to varying At. They get 
slightly worse, i.e. the fraction Vi is decreasing and the errors a{A.t) and o"(^) are increasing, 
for large time delays At, because the errors of the reconstructed lightcurves at the epoch of the 
interferometric observations increase with rising AtR due to a smaller number of error-reducing 
iteration steps needed for the reconstruction (cf. Figs. Q and Of course, the result also 
gets worse for very short values of At which are then becoming comparable to the timescale of 
variability and the observing interval AT. 
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Figure 10: The dependence of the results 
of the lightcurve reconstruction method 
on the time delay Ai. All other param- 
eters are kept fixed at the standard val- 
ues given in Tables ^ and ^ Plot (a) 
shows the variation of Pi (sohd line) , 1^2 
(dashed line) and Pa (dotted line). Plot 
(b) depicts the quantities M3 (solid line) , 
50 A/2 (dashed line) and Mi (dotted line). 
Diagrams (c) and (d) show la (solid 
line), 2cr (dashed line) and 3cr (dotted 
line) error estimates for the time delay 
and the magnification ratio. Note that 
the 3(7 lines are poorly determined, be- 
cause they are calculated from the AfR 
and /iR values for a few lightcurves only. 
(See the text for the definition of the 
quantities displayed in this figure.) 



15 




■ Figure 11: The dependence of the results 
'■ of the lightcurve reconstruction method 
on the magnification ratio /i. See the cap- 
j tion of Fig. ^ for the explanation of the 
1(J various lines in these plots. 



Magnification ratio ^u: FigureO verifies the statement about the magnification ratio which 
already has been given in Sect. |3.1| : The method completely fails for fj, 1 because then the 
errors due to the arbitrary initial guess do not decrease during the reconstruction process. In 
the diagrams the breakdown of the method for /i — > 1 is reflected by the steep decrease of Vi 
and the steep increase of the number Mj of minima appearing in the function. For /i ^ oo 
the vaules Vi, Mi and a{At) are remaining virtually constant. Note, however, that in this 
figure the relative accuracy of the flux density ratio measurements has been kept constant. For 
rising magnification ratio it will become increasingly difficult to measure the then also rising 
flux density ratios with the desired accuracy, and this will limit the applicability of the method 
for large 



Timescale of variability Tvar^ Figure 12, showing the dependence of the results on the 
timescale of variability, mainly reveals two effects. For very short timescales Tvar ^ 2 the errors 
introduced by interpolating the observed lightcurve between consecutive data points become 
large. This leads to an increase in the number of unwanted wrong minima and to a decrease 
of the fraction Vi- For slow variation, on the other hand, the number of minima roughly stays 
constant and Vi is only slightly decreasing, at least in the range of timescales plotted here. 
However, the error a{At) in the location of the true minimum is rising then. This means that 
the minima in the function are becoming broader and that the accuracy of the time delay 
determination is reduced. 




Figure 12: The dependence of the results 
of the lightcurve reconstruction method 
3& on the timescale of variability Tvar- 
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Figure 13: The dependence of the resuhs 
of the lightcurve reconstruction method 
on the variability u = as/{S\+2)- Here 
the noise parameter i] = ags /erg has been 
adapted such that the single-dish observ- 
0.00 O.OE 9.<H O.Oe i},Oe O.IO ing errors ass remain constant. 
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Variability v: A very important parameter is, of course, the relative variability of the source. 
In Fig. 13 not only u = as/{Si+2) has been varied, but also the noise parameter rj = ass fa's, 



which also includes the dispersion as of the combined lightcurve, has been adapted in such a 
way that the absolute errors ass of the single-dish lightcurve measurements remain constant. 
Naturally for u = no time delay determination is possible, but for increasing variability the 
results improve quite quickly, and Vi nearly reaches the ideal value of 1.0 for z/ > 0.04. The 
simulations have not been extended to variabilities larger than 10% because then the minima in 
the function are becoming very steep and narrow, which makes it more difficult to find them 
numerically. 



Number of interferometric observations A'^i: Figure 14 shows the evolution of the results 
when increasing the number of additional interferometric observations which are used as con- 
straints for the fit. From these diagrams it can be seen how the fraction of lightcurves for 
which the true minimum is the global minimum of the function rises and the number of 
minima and the errors in the location of the true minimum decrease. For eight fiux density 
ratio measurements, the fraction Vi has nearly reached 100% and typically there is just one 3a 
minimum remaining. 
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Figure 14: The dependence of the resuhs 
of the lightcurve reconstruction method 
on the number A'l of interferometric ob- 
servations. 
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Figure 15: The dependence of the resuhs 
of the lightcurve reconstruction method 
on the time interval ATi between the in- 
terferometric observations. 



Interval between interferometric observations ATj: For simplicity this study has been 
restricted to constant time intervals between the interferometric observations. These are per- 
formed towards the end of the observing period when the reconstruction is most accurate (cf. 
Fig. Figure 15 shows that the standard value ATj = 5 is almost ideal for the standard 
parameter set. The results get worse for small ATi which are comparable to the timescale of 
variability, because then consecutive flux density ratio measurements do not provide independent 
constraints. Therefore larger values of ATi are adequate for more slowly varying lightcurves. 
Very large values for ATj, however, are disfavourable, because then the first interferometric 
observations would have to be placed early on during the observing period when the errors of 
the reconstruction are still large. 
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Figure 16: The dependence of the results 
of the lightcurve reconstruction method 
on the observing period T. Plots (a)- 
(d) are of the same structure as in the 
previous figures. Diagrams (e) and (f) 
show the results for Vi and Mi for the 
same set of parameters as in (a) and (b). 
However, the interval considered as the 
range of potential time delay values has 
been reduced here from [0, 50] to [0, 30]. 
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Observing period T: In contrast to the autocorrelation method (cf. Fig. |2|a) the dependence 
of the results on the observing period T is rather indirect for the reconstruction method. The 
need for a sufficiently long observing period arises from the fact that the uncertainty due to 
the arbitrary initial guess gradually decreases during the reconstruction process. Eventually, 
the errors of the reconstructed lightcurves will be dominated by the observational errors of the 
single-dish measurements and will therefore roughly stay constant with time. This behaviour 



is reflected by the diagrams a to d of Fig. 16, For small T the results improve when extending 
the observation period, because then the interferometric measurements are still performed at 
epochs for which the reconstruction is affected by the initial errors. But since these have already 
dropped out at the respective epochs for T > 100, prolonging the observation period beyond 
that value does not improve the results any longer. 



Figure ICe and f show the effects of reducing the interval of potential time delay values, for 
which the reconstruction is done and for which the function is calculated, from [0, 50] to 
[0, 30]. Of course, this significantly improves the results, because then the wrong minima which 
are located in the interval ]30, 50] are excluded from the analysis. In addition, the worsening of 
the results for T < 100 is not so serious in this case, because these problems preferentially arise 
due to the large At^ values for which the errors of the reconstruction are decreasing slower, and 
just these At^ values have been excluded here. 

Observing interval AT: The effects of changing the observing interval AT for the single-dish 
measurements are conveniently demonstrated by inspecting the contours for an individual 
lightcurve realization. We have mentioned before that for calculating the reconstruction with 
AtR values which are not integer multiples of AT, the observed combined lightcurve has to be 
interpolated which in general causes additional errors for the reconstructed lightcurves. This 
leads to the effect that integer multiples of AT can be excluded as values for the time delay with 
a higher confidence than non-integer multiple values and that wrong minima are preferentially 



appearing at the latter ones. This is drastically demonstrated by the example shown in Fig. 17a. 
Here the observing interval has been changed to AT = 3 and a variety of wrong minima are 
showing up. Note, however, that the global minimum still corresponds to the true minimum. 
In practice the effect described here will be alleviated to some extent by irregular observing 



intervals. The example plotted in Fig. 17b shows that the adequate choice of the observing 



(a) 



(b) 




Figure 17; (a) This diagram shows the la and 3a con- 
tours of the function for an individual lightcurve 
realization with the standard parameter values, except 
that the single-dish observing interval has been changed 
to AT = 3. (b) For this example the observing interval 
is AT = 3 as well, but in addition to that the timescale 
of variability has been increased to Tvar ~ 15.9. 
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Figure 18: The dependence of the results 
of the lightcurve reconstruction method 
on the errors of the single-dish fiux den- 
sity measurements, characterized by the 
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interval sensitively depends on the timescale of variability. Despite using the same observing 
interval of AT = 3 as in Fig. 17a only one minimum is appearing here and the time delay can 
unambiguously be determined, because in this case the variability timescale Tvar = 15.9 is very 
much longer and therefore the interpolation errors are not so serious. 



Observing erros rj: From Fig. 18 it can be seen that fairly accurate single-dish flux density 
measurements are required for the method to work. For the standard parameter values the rms 
errors ass should not exceed 10% of the dispersion as of the lightcurve itself (ry = ass/o's = 0.1). 
The reason for the sensitivity to r] is, of course, that errors in the observed total lightcurve 
directly translate into errors of the reconstructed individual lightcurves and flux density ratios. 
In particular, it is desirable to determine the total flux density to a very good accuracy at the 
epochs of the interferometric observations, because the errors of these values contribute to the 
errors of the predicted flux density ratios to a large extent (cf. Eqs. (|l6|) and ([ITD). 



Flux density ratio erros a[m]/m: A crucial parameter is the accuracy of the additional 
constraints which are provided by the interferometric observations. Figure ^ depicts the de- 
pendence of the results on the relative error a[m]/m of the flux density ratio measurements. For 
the standard parameter values this error should not exceed ^ 1%. 
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Figure 19: The dependence of the results 
of the lightcurve reconstruction method 
on the relative error of the interferometric 
0.00 O.OB 0.04 O.Oe 0,08 0,10 flux density ratio measurements. 
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To complete the discussion of the parameter dependence we mention here that the effects of 
increasing the error (t[(5'i+2)] of the average combined flux density, which enters in Eq. (0), are 
negligible even for high values. Furthermore, changing the shape of the power spectrum which 
is used for generating the simulated lightcurves does not affect the results significantly. 



3.5 At determination for individual realizations 

Figure pO| a shows the contours for one example lightcurve generated with the standard pa- 
rameter values (Tables |^ and ^). In this case two minima appear with values low enough 
to be considered as candidates for the true minimum. Ideally one would like to include enough 
interferometric measurements so that only one sufficiently deep minimum is remaining. But 
here we have chosen this example in order to illustrate how we can determine probabilities for 
particular minima to be the true minimum in such cases and how error ranges for the time delay 
and the magnification ratio can be obtained. The former will not be possible without making 
use of the simulations, but concerning the latter point we can again apply statistics. 
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Figure 20; (a) The 3a contours (as defined in Sect. 3.3) 
for one iigfitcurve realization witli tlie standard param- 
eter values. The cross marks the global minimum. Plot 

(b) shows the 3cr, 2a and la confidence regions for 
At and /i under the assumption that the minimum at 
AtR — 19.9 is indeed the true minimum. Accordingly 

(c) shows the confidence regions for assuming that the 
minimum at AtR — 36.3 is the true minimum. The lim- 
iting values for these confidence regions were derived 
from the x^-distribution with 2 degrees of freedom for 
Ax^ = Xexact - XtTuc (sec Fig. |lj). Thus, diagrams (b) 
and (c) show the contours of x^ = Xtruc + with 
Ax^ taken from the vertical lines in Fig. ^ and Xtruc 
taken as the value at the respective minimum. 



Fi gure |21| verifies that the distribution of the difference — Xexact — Xtrue between the 
values at the exact parameter values {At, /i) and at the true minimum is following a x^- 
distribution with 2 degrees of freedom, as it is theoretically expected for two fit parameters {At^ 
and /xr). Assuming that a particular minimum is indeed the true minimum, we can use the 
limiting Ax^ values indicated in Fig. |2T| to obtain confidence regions for the time delay At and 



the magnification ratio /x. This has been done in plots c and d of Fig. 2C for the minima at 
AtR = 19.9 and At^ = 36.3, respectively. Note that for determining confidence intervals for the 
time delay At alone, irrespective of the magnification ratio, the x^-distribution with 1 degree of 
freedom has to be used for Ax^. 

In order to quantify the probability for the respective minima to be the true minimum we 
distinguish the following two cases: 

A: The minimum at At^ = 19.9 is the true minimum. Then the minimum at At^ = 36.3 
necessarily has to be the lowest wrong minimum. 
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Figure 21: The histogram shows the distribution 
of Ax^ = xLact - Xtrue determined from the sim- 
ulated lightcurves (with standard parameters) that 
have already been used for generating Figs. ^ and 0. 
This distribution is consistent with the theoretical x " 
distribution with 2 degrees of freedom which is plotted 
in the diagram as well. According to this distribution 
the dotted vertical lines indicate the limiting Ax'^ val- 
ues for a given exclusion confidence. 



B: The minimum at Atp;, = 36.3 is the true minimum. This means that the minimum at 
AtR = 19.9 is the lowest wrong minimum. 

Here we are neglecting the possibility that both of these minima are wrong, because this would 
imply that the true minimum is at a \^ value beyond the 2>a limit and thus very far in the 
exponentially decreasing tail of the probability density distribution (cf. Fig. In Fig. [2^ we 
again show the probability density distribution Ptvueix^) the true minimum, represented by 
the theoretical ^^-distribution with 3 degrees of freedom, and the probability density distribution 
Pwrong(x^) for the lowest wrong minimum as determined from the simulations. In the following 
we denote the values at the A^r = 19.9 minimum with Xa ^-^d those at the At^ = 36.3 
minimum with Xb- the figure these values are marked with vertical lines, the solid line for 
1 and the dashed line for Xb ^ 6. In case A, i.e. At^ = 19.9 being the true minimum, the 



X, 



value for Xa is drawn from the distribution Ptrue(x^) ^-nd therefore Xb is drawn from Pwrong(x^)) 
whereas for case B it is vice versa. In Fig. the values Xwrong ^t the lowest wrong minimum 
have been plotted against the value Xtrue ^t the true minimum for the individual realizations of 
the lightcurve simulations. Since this scatter plot shows no correlation between these x^ values, 
we can treat Ptrue(x^) and Pwrong(x^) as independent distributions. This allows us to assign 
probabilities to the cases A and B: 



Ptr 



3(Xa)Pwrong(Xb) 



J^true (Xa) Pwrong (Xb ) + Ptrue (Xfa ) Pwrong (X 



91% 



Ptruc(xI)Pwrong(Xb) + Ptruc (Xb) Pwrong (xl) 

Concerning the determination of the time delay and the magnification ratio we can therefore 
make the following statements. With a probability of 91% the minimum at A^r = 19.9 is the 



^Note that the 3a limit which has been used for the contours in Fig. |2C|a and which is indicated by the vertical 
line in Fig. ^ is derived from the x^-distribution with Ni — 5 degrees of freedom and therefore is conservative for 
the purpose of excluding wrong minima. To do th at th e distribution with A^i — 2 = 3 degrees of freedom would 
be adequate. See the discussion at the end of Sect. 3.3. 

Figure 22: (a) This plot again shows the 
distributions of the x^ values at the true 
minimum, represented by the theoretical 
curve (solid line), and at the lowest wrong 
minimum (dashed line), determined from 
the simulations (cf. Fig. 0). The thick 
vertical lines indicate the x values at the 
two minima of the function shown in 
Fig. §o| (solid line for AtR = 19.9, dashed 
line for AtR — 36.3). (b) In this scatter 
plot the x^ values at the lowest wrong 
minimum are plotted against those at the 
true minimum for each realization of the 
simulated lightcurves. 
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true minimum. Then confidence regions for tlie time delay At and tlie magnification ratio /x are 
given by Fig. |20|b. Tlie probability for Atpj, = 36.3 being the true minimum is 9%. In this case 



the confidence regions for At and /i are given by Fig. 203. In view of the much higher probability 
of case A, which is indeed the true one, this result could also be phrased as: The error level for 
restricting the time delay to be At = 19.9 with confidence regions as shown in Fig. 2Cb is 9%. 

At this point a few cautionary remarks should be added. The probabilities derived here rely 
on the distribution of the wrong minima which cannot be described by statistics and therefore 
has to be determined from the simulations. Hence, in applications to real observations it is 
important to generate a set of simulated lightcurves which satisfactorily reflects the variability 
characteristics of the observed lightcurve. Determining the confidence regions in Fig. 20 b and c, 
however, solely relies on statistics and does not require any simulations. Of course, for any 

fit it is essential to use error estimates which are really reflecting the statistical uncertainties 
of the measurements and not, e.g., the observer's predilection for giving conservative limits. 

To conclude this section we demonstrate how the results can be improved by additionally 
taking into account the information contained in the autocorrelation function. In Fig. |23| a the 
value C(r)|„rong of the autocorrelation function for the lag r corresponding to the time delay 
value AtR of the lowest wrong minimum has been plotted against the value at that minimum 
for the individual lightcurve realizations. The scatter plot indicates that there is no correlation 
between these two quantities and that on average the value C(t)| wrong is zero, i.e. there is no 
enhanced autocorrelation for r values corresponding to the At^ values of wrong minima. This 
implies that the lightcurve reconstruction method and the autocorrelation method are completely 
independent techniques for determining the time delay. Figure 23b depicts the distributions of 
the autocorrelation values for lags r corresponding to the At^ value of the true minimum (solid 
line) and of the lowest wrong minimum (dashed line). The latter distribution is centered around 
C{t) ~ 0, whereas the former is shifted to higher autocorrelation values. The reason for this 
is, of course, that the true minimum is roughly located at the correct time delay At and so 
the effect discussed in Sect. § leads to an enhanced autocorrelation. From the large overlap 
of these distributions and from the small separation of the vertical lines in this plot, which in 
analogy to Fig. 22 a denote the values for the two minima of the example case, it is clear that 
the autocorrelation function cannot improve the results drastically. Nevertheless, including this 
additional information in a similar analysis as above leads to probability values of Pa ~ 95% 
and Vb ~ 5% and thus the error level for tying oneself down to At = 19.9 as the value for the 
time delay could be reduced from 9% to 5%. 



(a) 1.0 




Figure 23; (a) In this diagram the val- 
ues of the autocorrelation function for 
lags T corresponding to the AtR value of 
the lowest wrong minimum are plotted 
against the value at that minimum, 
(b) This plot shows the distribution of 
the C(r) values for r corresponding to 
the true minimum (solid line) and to the 
lowest wrong minimum (dashed line). In 
analogy to Fig^ 2^ a the values at the min- 
1.0 ima of Fig. ^ have been indicated by 
thick vertical lines. 
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4 Discussion and conclusions 



In the first part of tliis paper we have shown that an analysis of the autocorrelation function 
of the total lightcurve of gravitational lens systems can in principle be used to determine its 
time delay. In practice, however, a fairly long observing time is required to make sure that the 
measured autocorrelation function reliably reflects the lightcurve's underlying statistical prop- 
erties. In addition, the characteristic time scale for the source variability has to be considerably 
shorter then the time delay, and the magnification ratio should be close to unity in order for the 
autocorrelation function to develop distinct time delay features. 

It is interesting to note that it has been proposed by Press (1996) that in certain cicumstances 
time autocorrelation effects on unresolved, time- varying objects could be used to establish other- 
wise undetectable (micro-) lens systems. In view of the results shown in Sect. ^ a detection of the 
time delay effect via an analysis of the autocorrelation function seems to require a delicate fine 
tuning between the involved time scales, i.e. the observing interval, the time scale of variability, 
the time delay, the observing period and possibly the "life time" of the lensing configuration. 
However, using different approaches might be more promising, especially for non-gaussian ran- 
dom processes. Consider, as an extreme example, an intrinsic source lightcurve consisting of a 
characteristic series of sharp peaks. In such a "burst-scenario" a time delay could be easily de- 
termined from the repetition of the features which is then showing up in the observed lightcurve 
of the lens system. 

In Sect. 1^ we introduced the "lightcurve reconstruction method". We demonstrated that 
it is possible to reconstruct the individual lightcurves of a gravitational lens system from the 
observed total lightcurve by assuming values for the time delay and the magnification ratio. 
The reconstructed lightcurves are most accurate at the end of the observing period, because er- 
rors introduced by an arbitrary initial guess are gradually decreasing during the reconstruction 
process. However, this is not true for a magnification ratio close to unity and the method is 
therefore not applicable in such cases. In general it will be impossible without additional infor- 
mation to single out the true time delay from all possible reconstructions with potential values 
for the time delay and the magnification ratio. (A counterexample for which this might indeed 
be feasible by assuming some reasonable shape for the correctly reconstructed lightcurves is the 
highly non-gaussian burst-scenario outlined above.) However, interferometric measurements of 
the flux density ratio (at the same frequency as the total flux density monitoring!) provide 
an excellent means for checking the consistency of the various reconstructions and eventually 
determining the time delay unambiguously. 

With simulations we investigated the dependence of the method on various parameters. 
Naturally the variability of the lightcurve turns out to be an essential parameter - as it is 
the case for any time delay determination technique. In addition, a fairly high accuracy for 
the total flux density monitoring observations as well as for the additional interferometric flux 
density ratio measurements is desirable. In contrast to the autocorrelation method the lightcurve 
reconstruction also works for rather slowly varying lightcurves, although increasing the time 
scale of variability leads to larger confidence intervals for the time delay. (That lightcurves 
with fairly slow variability still contain information about the time delay can also be seen from 
van Ommen et al.'s (1995) analysis of the lightcurve of PKS1830-211, a lens system with two 
compact images and an Einstein ring. They use simple linear fits to the slowly varying total 
lightcurve and combine these with interferometric flux density ratio measurements to get an 
estimate for the time delay. However, the confidence intervals for this value are still quite large.) 
Of course, the simulations that we have performed just show the general trends and cannot 
cover the parameter space completely. For a given lens system and observation programme 
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some unfavourable parameter values could be compensated by another parameter being much 
more adequate and thus leading to equivalently good results. 

In principle the method described here could be used for radio as well as for optical lightcurves. 
However, there are several arguments against optical wavelengths. If the time scale of variabil- 
ity is short, a small observing interval is essential to sample the lightcurve accurately and a 
larger gap due to bad weather would have drastic effects on the quality of the reconstruction. 
In addition, microlensing effects which cannot properly be acounted for in this method can be 
important in the optical. However, if the microlensing time scale is considerably larger than the 
time delay and the observing period, this will only result in a change of the effective magnifi- 
cation ratio and the method can still be applied. Finally, a very important reason to perform 
such observations at radio frequencies is that the most promising lens systems for determining 
the Hubble constant consist of double images of radio-loud sources and an associated ring of 
extended emission. From these Einstein rings the system B0218-I-357 mentioned in the introduc- 
tion probably is the best candidate at the moment. The time delay for this system is expected to 
be roughly between 8 and 20 days which is reflecting our current knowledge (or rather ignorance) 
of the Hubble constant and uncertainties in the lens model. (A detailed model for this system 
which takes into account the information provided by the Einstein ring still has to be made.) 
The preliminary value of At = 12 it 3 for the time delay obtained by Corbett et al. (1996) 
is consistent with this range. Reliable estimates for the time scale and the amplitude of the 
source variability are required in order to design an optimal observing strategy and to specify 
the number of interferometric observations needed for a definite determination of the time delay 
by using the lightcurve reconstruction method. However, combining a single-dish monitoring 
programme of two or three months with five to ten additional interferometric constraints should 
be adequate to determine the time delay reliably. 

For the system PKS1830-211 the discovery of molecular absorption lines by Wiklind & 
Combes (1996) opens an interesting new possibility to determine the time delay from unre- 
solved single-dish monitoring observations only. Apparently only one of the two source images 
is covered by a molecular cloud, and this allows one to determine the flux density ratio from 
molecular spectroscopy. Combining this with the total flux density it is possible to calculate 
the individual lightcurves and apply a cross-correlation analysis in order to obtain the time de- 
lay. However, the lightcurve reconstruction method could be used in this case as well, with the 
advantage of a much cleaner error treatment, because it uses the direct observables total flux 
density and flux density ratio. 
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